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General relativistic corrections to the expansion rate of the Universe arise when the Einstein 
equations are averaged over a spatial volume in a locally inhomogeneous cosmology. It has been 
suggested that they may contribute to the observed cosmic acceleration. In this paper, we propose 
a new scheme that utilizes numerical simulations to make a realistic estimate of the magnitude of 
these corrections for general inhomogeneities in (3+1) spacetime. We then quantitatively calculate 
, the volume averaged expansion rate using N-body large-scale structure simulations and compare it 

. with the expansion rate in a standard FRW cosmology. We find that in the weak gravitational field 

limit, the converged corrections are slightly larger than the previous claimed 10 -5 level, but not 
large enough nor even of the correct sign to drive the current cosmic acceleration. Nevertheless, 
the question of whether the cumulative effect can significantly change the expansion history of the 
Universe needs to be further investigated with strong-field relativity. 
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I. INTRODUCTION 
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JL ' One of the most puzzling questions in modern cosmology is the nature and origin of the dark energy that is 
; i . responsible for the present cosmic acceleration. Over the last decade, evidence has been accumulating from the type 
' la supernova luminosity distance-redshift relation [1-3] and other observations [4-8] indicating that the Universe 
C3 '. is accelerating. In the standard Friedmann-Robertson- Walker (FRW) cosmology model, this can be explained by 
introducing a dark energy term in the Friedmann equations. Observations indicate that the dark energy comprises 
[ more than 70% of the mass-energy in the Universe. Numerous explanations have been proposed for the origin of the 
■ dark energy. The simplest explanation is that of a mass-energy that violates the strong energy condition either in 
the form of a cosmological constant [1-3] or a quintessence [9, 10]. However, this solution has unnatural hne tuning 
, and coincidence problems. Alternately, one could assume that GR is not the complete theory on cosmological scales 
(such as in the Dvali-Gabadadze-Porrati (DGP) [11] or f(R) [12, 13] models), however, no corroborating evidence for 
deviations from standard GR has yet been found. An alternative explanation of interest to the present work is that 
£NJ ■ corrections to a FRW cosmology, due to the presence of local inhomogeneities, may introduce dark energy like terms 

' in the equation of cosmic expansion. 
On . While the Universe appears homogeneous and isotropic on large cosmological scales, this is not the case on smaller 
scales. Local inhomogencity and structure always exist. It has been argued [14-16], therefore, that one should 
study the observational data in the context of a realistic lumpy universe first, instead of assuming that the FRW 
model is correct, to fit cosmological parameters. Hence, it is imperative to clarify how the inhomogeneities affect the 
interpretation of the- observational data. 

Even before the discovery of the cosmic acceleration, it was argued [17] that when light propagates only through 
a nearly empty intergalactic medium, there is a dimming effect compared to the FRW model. This may affect the 
interpretation of the supernova luminosity distance-redshift data. Recently, one simple class of the inhomogeneous 
models, the Lemaitre-Tolman-Bondi (LTB) model, has been extensively studied [18-21]. In this model, it is assumed 
that we are living near the center of a spherical underdense region and the Universe may have one or many such 
regions. The supernova data can then be fit without introducing a dark energy term. Although we may indeed 
reside in an underdense region, the special symmetry in this model is not consistent with current observations of the 
large-scale structure of the Universe. Moreover, none of the realistic models of this class can fit all of the observed 
cosmological constrains [19, 20, 22, 23]. 
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While more general inhomogeneous cosmological models have been investigated using various methods [24-27] , the 
domain averaging procedure proposed by Buchcrt [28-31] has been of most interest recently. By using perturbation 
methods in this domain averaging procedure in a general synchronous gauge, Kolb et. al [32, 33] claimed that an 
effective negative pressure term can arise from averaging local fluctuations. However, it has been strongly argued 
[34-39] that the effect from this procedure is either not a physical observable or not large enough to drive the current 
cosmic acceleration. Moreover, such perturbativc analysis is not easily amenable to the nonlinear evolution of the 
large-scale structure. The lack of a realistic model and the limitation of the perturbation methods make it difficult to 
calculate the proposed effect accurately. To alleviate this problem, therefore, we describe here the first step toward 
the development of a relativistic numerical scheme to explicitly calculate the proposed effect of domain averaging. 

The paper is organized as follows: We first develop a theoretical scheme that enables us to calculate the magnitude 
of the proposed correction terms quantitatively using an N-body large-scale structure simulation as described in Sec. 
II. The details of the numerical simulations are discussed in Sec. III. The results are presented and discussed in Sec. 
IV followed by a summary in Sec. V. The main purpose of this paper is to present the formalism and to make initial 
numerical investigations of the effects from this domain averaging procedure on the expansion rate of the Universe. 
This study will shed light on whether this is a viable approach to explain the nature of the dark energy 



II. GENERAL RELATIVISTIC CORRECTIONS FROM DOMAIN AVERAGING 



A. Domain averaged expansion rate in an inhomogeneous cosmology 

Here, we will roughly follow the domain averaging procedure proposed by Buchert [28-31]. Numerical calculations 
in general relativity are best formulated in the Arnowitt-Deser-Misner (ADM) formalism [40-42]. Hence, we start 
with the general (3+1) ADM metric: 

ds 2 = -(a 2 - Pifi^dt 2 + 2fcdx i dt + -/ lj dx l dx j , (1) 

where a is the lapse function denoting the lapse of proper time, pi is the shift vector describing the shift of coordinates, 
respectively, from one time slice to the next and 7^ is the spatial three-metric. The proper volume of an arbitrary 
domain D in this scheme can then be defined as: 

V D = f jd 3 x , (2) 



where 7 = y/ det(jij ) , and detfoij) is the determinant of 7^. We can then define the average of an arbitrary scalar 
field ip(x,t) on the domain D as: 

(l/>(x,t)) D = ^r I ^(x,t)-yd S X . (3) 



The time derivative of this domain average is then: 

d(tp)p d 1 



Now, using the fact that for any invertible matrix A 



= -^> D + (4 + F / Tp7d 3 x . (4) 



d -^m = det(A)(A-i 



detiA^A- 1 )^ , (5) 



and the fact that 7^- is symmetric. The derivatives of 7^ can be written as: 
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We next choose the Eulerian gauge (shift vector = 0)* to reduce the ADM metric to: 

ds 2 = -a 2 (x, t)dt 2 + jij(x, t)dx l dx 3 . (7) 

In this gauge, the extrinsic curvature Kij , which can be interpreted as the rate of change of the spatial metric 7^ 
along the normal vector, can be simply expressed as: 

Ku = ~ia, K = ^K ij = ~y ii ^ . (8) 

la la 
By inserting the expression for |7y7 y from Eq. (8) into Eq. (6), we have: 

7 = —aKj . (9) 
The time evolution of the domain average in Eq. (4), can then be expressed as: 

Vd -Wo + Wo-^- f aKvixPx . (10) 



dt V D v Vd j 

The new dimensionless scale factor a D for the domain D, and the new Hubble parameter H D for this domain, can 
now be defined as: 

^V /3 , H D = ^ = l^, (11) 
V Do J a D 3 Vd 

where V D is the proper volume defined in Eq. (2) and V Do denotes the domain volume at the present time. H D can 
be further expressed as: 

where we introduce the trace of the expansion tensor O*- as Q = Q\ = —aK. H D is now the appropriate physical 
quantity to describe the overall expansion rate of the inhomogeneous domain D. In what follows, we calculate this 
quantity explicitly using numerical large-scale structure simulations as realistic representations of the evolution of the 
lumpy universe within the domain. 



B. Time evolution of the expansion rate and effective pressure 

We now derive the time evolution equation of the domain averaged expansion rate. Substituting Eq. (12) into Eq. 
(10), we obtain an important commutation rule: 

(Q)»(4>)n , (13) 

-<9} 2 D ■ (14) 
In the ADM formalism, the time evolution of the extrinsic curvature scalar K can be written as [40-42] : 

K= -D i D i a + a( 3 R + K 2 ) + 4nGa(S -3p H ) , (15) 

where D t is the covariant derivative operator in the three-space which reduces to an ordinary gradient operator for a 
scalar such as a. The quantity 3 R is the Ricci scalar for the three-metric 7^. S = 3P + ph(W 2 — 1) is the trace of the 



which for %jj = 0, gives: 

m ° -<e) D = (e 2 ) D 



m 



* The purpose of this simplification and application of the conformally flat condition in Sec. IIB is solely to make easy use of current 
numerical large-scale structure simulation codes. The generalization of the expressions in Sec. IIA and IIB in the general ADM formalism 
is straightforward. 
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spatial stress, and pn = phW 2 — P is the Hamiltonian density. In the expressions above, p is the rest mass-energy 
density, P is the pressure, W = vT + uiu 1 is a generalized Lorentz factor, m denotes spatial components of the four 
velocity, h = l + e + P/pis the specific enthalpy, and e is the internal energy per unit mass. 

The Hamiltonian constraint derives from the "00" component of the Einstein equation. In the ADM formalism, it 
can be written as: 

3 R + K 2 - K,,K ,! = 16nGp H ■ (16) 
In the slow-motion, low-temperature limit, i.e. W = 1, e = 0, pn = p, S — 3pH — 3(P — p), Eq. (15) then becomes: 

K = -DtD'a + a( 3 R + K 2 ) + l2nGa{P - p) . (17) 

Henceforth, we denote 3 i? as R. 

For the present application, there is no significant rotation, shear or gravity waves. Hence, we can adopt a confor- 
mally flat condition [42] to further simplify the ADM metric to the form: 

ds 2 = -a 2 {x,t)dt 2 + a 2 FRW {t)(l) 4 {x,t)6 i:j dx l dx j , (18) 

where a FRW (t) is the scale factor in the FRW model limit which is generally different from a D in Eq. (11), <p(x,t) 
is the conformal factor denoting the local deviations from a homogeneous and isotropic curvature in the three-space. 
Now, the conformally flat three-metric 7^ takes the form: 

7ii=4 w (^ 4 (a;,^i ■ (19) 
In this formulation, Kij has diagonal elements only and KijK 1 ^ = —^K 2 , thus we have: 

la a a FRW q> 

K = _l { o^ L + 2 k (20) 

a a FRW 

The Hamiltonian constraint can then be reduced to: 

R + \k 2 = WnGp . (21) 

Using Eq. (21), the K equation [Eq. (17)] can be rewritten as: 

K = —D i D i a + \aK 2 + 4nGa{p + 3P) . (22) 
Now employing the relation aK = —0 — dK 7 we can express the domain averaged K equation as: 

(ak) D = -(e) D -(aK) D = -(aD l D l a) D + ^(e 2 ) D +4irG(a 2 (p + 3P)) D . (23) 



Using Eq. (14), we can then rewrite Eq. (23) as: 



(6 2 ) D - (0)| - ^£ = ( (aK) D - (aDiD'a)^ + ho 2 ) D + AvG(a 2 {p + 3P)) D . 
at 3 

The d(Q) D /dt term can be rewritten as: 

^ = 3 ^^ = 3 /a £ _/a £ \ 2 \ =3 a £ _l % 
at at \a D \a D J J a D 3 

Plugging this expression into Eq. (24), we have: 

3 ^ + S D ( a )+4irG(a 2 (p + 3P)) D = l((e 2 ) D -(e)l) . (26) 

(In O 
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Here, wc denote S D (a) = (aK) D — (aDiD l a) D because it behaves like a source term. The new term | ((6 2 )u — (©)?>) 
is the kinematic backreaction term Q D originally proposed by Buchert [28, 29]. 

Now, we can write the modified first Fricdmann equation by introducing an effective energy density p e ff, 

H ° = {t) 2 = r Gp ° s ' with ^ = 2ib (e)2 ° ■ (27) 

Similarly, from Eq. (26), the modified second Friedmann equation can be written by introducing an effective pressure 

PeS, 

^ = -^(Pe ff + 3P off ) , with P cff = I(a 2 (p + 3P)) D + I ^5 I3 (a) + ^<e) 2 D - T ^(e 2 ) D . (28) 

Here, = — aK = 3(H FRW + 2<p/(f>) and 2 = a 2 K 2 = 9(H FRW + 2<p/(j)) 2 , where H FRW is the Hubble parameter 
defined by the unperturbed Fricdmann equation, H FRW = (a FRW / a FRW f = l>KGp FRW and 6 -> 3H F RW in the FRW 
limit. We will use the correction from p FRW to p c s to represent the correction from H FRW to H D in Sec. IV. 



C. Evaluating the correction terms with large-scale structure simulations 

In order to connect the correction terms that we derived in Sec. II A and IIB with presently available large-scale 
structure simulation codes, we begin with the conformal Newtonian gauge [43, 44], 

ds 2 = a 2 FRW (r 1 )[-{l + 2<5>)dr l 2 + {\-2>$>)dx % dx l ] . (29) 

This describes a restricted class of general gauge-invariant cosmological perturbation theories [45-48] . We then identify 
this as the weak-field limit of our conformally flat metric [Eq. (18)]. Here, $ is the peculiar gravitational potential, r\ 
is the conformal time, drj = dt/a FRW (t). Using the fact that for any given rj, there is a corresponding t, we can make 
a mapping that a 2 — > 1 + 2$ and </> 4 — > 1 — 2$ to express the metric coefficients, a and </>, in terms of $. 

By definition, the domain averaged expansion rate H D is a constant within a local domain at a given redshift. 
Utilizing this, and the expression <j)/<j) = — $/2(l — 2$), H D can be expressed as: 

H D = \{-aK) D = (H FRW + 2^) D = H FRW + 2{^) D = H FRW - {-^rr)o ■ (30) 
6 (f> (f> 1 — 2$ 

The determinant of 7^ in Eq. (19) now becomes, 7 = a FRW (t)(f> 6 (x, t) = a FRW (l — 2$)^, and by using the definition 
of the domain average in Eqs. (2) and (3), H D can be further written as: 

LA^d 3 x f 1(1-2$)?^ , 

77 TT J P 1 — ' JT J D V ' 

n D — n FRW f . — n FRW 3 . {01 

J D jd 3 x f n (l-2$)2d 3 x 

Now, H D is only dependent on H FRW and $. Both quantities can be easily extracted from large-scale structure 
simulations. The effective energy density p e ff can be calculated as: 

_3Hl_ 3 u2 _3_ f <j>(l-2<#rf 3 x 3 f J D Ml-m^d 3 x \ 2 

PcS ~ SttG ~ ^G Hfrw AkG Hfrw J d (1 _ 2$) I d?x SttG { J d (1- 2$) I d 3 x ) " [ ' 

The correction to the energy density, (p c g — p FRW ) / ' p FRW i can then be expressed as: 



Peg - Pjmw = 2 f D *(l-2*)i<Px _J_ ( J D Ml - 2^ d 3 x 

p FRW H FRW J D (l-2$)§d 3 x H 2 mv I JJl - 2<S>)I d 3 x 



2 



(33) 



The expression for the effective pressure P e ff in Eq. (28) is somewhat cumbersome. Hence, we will analyze it term 
by term. Since for the current application, we only need to deal with a universe dominated by nonrelativistic matter, 
we can take the matter pressure P to be negligible. The first term in the expression of P c g in Eq. (28) can then be 
written as: 

1, „, S1 if j o(l + 2$)(l-2$)id 3 x / N 

— (a(p + 3P)) D = — ^-2— !± — 1 . 34 

3 V KF "° 3 J (l-2$)id 3 x V ' 
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If we define the density fluctuation as 6p = p — Pfhw 

and use the fact p FRW — 3H frw /St:G, we have: 

1,2, , 1 2 1 a ( $(l-2$)td 3 a; lf n Wl-4f 2 )(l-2$)i ( | 3 3: 

-(a 2 (p + 3P)) g = + -7ZrHl R J D ' + - Jp ^ ~ • (35) 

3 87rG 47rG J d (1 - 2$)5 d 3 x 3 J D (1 - 2$) 5 d 3 x 

For the second term, we have a = 6/(1 + 2$)i, and D^'a = -'■ ! l),l),<\ = 3(a| w (l - 2$))- 1 V 2 (l + 2$)*. The 
entire second term can then be written as: 

1 1 / $ $ i V 2 (l + 2$)2 \ 

~ = -— (——(H FRW --—)) D + ((l + 2^ 



12ttG uk ' 4ttG^1 + 2$ v 1-2*"" ^ v ' a% w (l - 2$) 

= -4^ /d( i - 2*)* cPx U " " 2<&)§ ^ 

+ y"(l-4$ 2 )^a; 2 w V 2 (l + 2$)^d 3 .Tj . (36) 

Using the fact that $ < 1 (or <£>/c 2 -C 1, if we explicitly denote the value of "c") in the weak- field limit and the 
cosmic Poisson equation: 

V 2 $ = ATTGa 2 FRW 5 P , (37) 
we deduce that a~ 2 w V 2 (l + 2$) 2 = AnGSp. So the second term can be further written as: 

1 q . , 1 „ ^(1-2$)* (fir l /^(l^)^ JJp(l-A^)U^ 

dn(a) - -— — M FRW — — , ,„ r ~ A 77 m nli 3 z, r— ^.3 ,„ • (38J 



12ttG " v ' 4ttG / D (l-2$)ld 3 x 4ttG fjl - 2<P)U 3 x / D (l-2$)iA 

The third term is simply the p c g. The fourth term can be written as: 

1 . (e2s _ 1 L(^w- T ^ 5 ) 2 (l-2$)l^ 



18ttG n 2ttG / (1-2$) id 3 ! 



1 i J fl ^-^)^fe _ i ^ fljjjT ^ 



2ttG * kw ttG *~" J d (1-2$)jA 2nGJ o (l-2$)i(Px 

The effective pressure P e ff can now be summarized as: 

1 - / p $(l-2$)fd 3 x 1 6^(1 - 2^)^(1 + 2c&)-i 

nS - ^G Ufrw f o (l-2*)$<Px + ttG Ufrw / D (l-2$)ld 3 x 
1 / D 5/9(1 - 4$ 2 )(1 - 2$)5 d 3 a; J d 8p(l - 4$ 2 ) 2 d 3 x 



3 / D (l-2$)ld 3 x ^(1-2$)!^ 

1 J £) $ 2 (l + 6$)(l + 2$)- 1 (l-2$)-5d 3 .T 3 f J D $(1-2$) 3d 3 



.i' 
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4ttG J D (l-2$)td 3 a; 8ttG ^ Jjl - 2$)f d 3 x / 

The equation of state parameter of the effective dark energy like term can be defined as w c g = P e s/(p e s — Pfhw)- 



(40) 



III. DETAILS OF THE NUMERICAL SIMULATION 



In the present work, our goal is to estimate the magnitude of the deduced correction terms. To achieve this, we 
do a straightforward large-scale structure simulation in a standard FRW cosmology and post-process the simulation 
data to evaluate the correction terms at each given redshift. In subsequent works, we will evolve the simulation with 
the modified equations of motion in real time and evaluate the resulting cumulative effect. 

The code we have adopted for the present numerical simulation is the N-body SPH code GADGET which was 
originally developed by Springcl ct al. [49]. The most current publicly available version, Gadget-2 [50], is used for all 
of the simulations described in this paper. 
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We set up the initial condition for the simulation as follows: First, the initial linear matter power spectrum for one 
specific set of cosmological parameters is generated by standard CMB codes such as CMBFAST [51]; We then use 
the associated power spectrum to generate a Gaussian random field with the Zel'dovich approximation [52] utilizing 
such generator codes as the Grafic [53] packages and the IC [54] package; Finally, the file is converted into the 
Gadget-2 format to start the simulation within a periodic comoving box. Some authors [31, 55-57] have argued that 
the application of a periodic boundary condition in numerical simulations is equivalent to forcing a standard FRW 
cosmology and thus eliminates any proposed corrections. However, for the tree algorithm in the Gadget-2 simulation 
[49, 50, 58, 59] and the size of our domain, this is not necessarily the case as discussed below. 

In order to calculate the new domain averaged terms, p e ff, -Peff and u> e ff , as derived in Sec. II, we utilize the Gadget-2 
output of particle masses and positions at each time slice. For the domain averaging procedure, we first divide the 
whole domain, i.e. the entire simulation box, into a fine spatial grid of size L^. We then use a Cloud-in-Cell (CIC) [60] 
method to assign the matter density p to each zone within the grid. Finally, we use the cosmic Poisson equation [Eq. 
(37)] with the method of Successive Over Relaxation (SOR) [61] to calculate the peculiar gravitational potential $ for 
each zone. We note that the application of a periodic boundary condition in the simulation implies that the solution 
to Eq. (37) for the whole grid is only unique up to an arbitrary constant, this is adequate for evolving the equation 
of motion because they only involve the gradient of $. However, this periodic boundary condition approach does not 
determine the value of $ uniquely as is needed to calculate the correction terms. Hence, we use a fixed comoving 
boundary assuming that outside of the simulation box, the matter density is exactly p FRW . In this way, $ is uniquely 
determined by the matter density distribution within the simulation only. Once p and <J> are obtained in this way, the 
domain averaged quantities defined in Sec. II can be calculated. The physical implication of this procedure is that 
by reducing the smoothing length from the Hubble scale to the resolution limit of the simulation, we can effectively 
evaluate the effect of the local inhomogencities on the global cosmic expansion rate, which is always neglected in the 
standard FRW cosmology. This is the procedure that was proposed by Ellis [14, 16] and by Ellis and Stoeger [15]. 
The accuracy of this approach is limited, however, by the number of particles in the simulation and the resolution of 
the spatial grid. In the next section, we examine the dependence of the effect on the size, resolution and number of 
particles in the simulation. We show that a converged result can be obtained. 



IV. RESULTS 



As an illustration, Fig. 1 shows the matter power spectrum for three different cosmologies as labeled from our 
simulations at redshift z=0. On small scales, the simulations involve substantial nonlinear growth of structures that 
a linear theory cannot predict accurately Fig. 2 shows the growth of the peculiar gravitational potential $ during 
the structure formation epoch at three different redshifts in a flat, Q m — 1 cosmology simulation. In all of these 
simulations, the upper limit of the absolute value of the metric perturbations, $/c 2 , is about 10 -3 and the upper 
limit of the magnitude of the peculiar velocities is about 10 4 km/s (v/c < 0.03). Hence, the simulations indeed stay 
within the weak-field, slow-motion regime as assumed in the derivation in Sec. II. 



a =0.3fl A =0.7 -- -- 

fi_=0.3 a.=o.o 

Q m =1.0ii A =0.0 

£i m =0.3 Q A =0.7 linear 

£i„ =0.3 Q A =0.0 linear 

£i„=1.0 ii.=0.0 linear 
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FIG. 1: Matter power spectrum for three different cosmologies from both a linear theory and simulations. 



Next, we study the proposed correction terms in detail using simulations in a flat, Q m = 1, matter-dominated 
cosmology. Table I lists the various box sizes and numbers of particles we used for this special case. 

Since the proposed correction terms arise from the domain averaging procedure, we wish to investigate their depen- 
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FIG. 2: The growth of the peculiar gravitational potential $ from redshift z = 9to2 = 0ona comoving spatial slice in the 
X-Y plane in a three dimensional simulation. The value of $/c 2 is drawn along the Z axis as indicated at the bottom of the 
figure. 



TABLE I: List of the parameters of the simulations. 



Parameter 


I 


II 


III 


IV 


V 


VI 


VII 


VIII 


IX 


Box size (Mpc/h) L 


100 


100 


200 


200 


400 


400 


800 


800 


1600 


Number of particles N p 


128 3 


256 3 


128 3 


256 3 


128 3 


256 3 


256 3 


512 3 


256 3 



dence on the smoothing length Ld, which is the individual zone size of the grid.^ Figs. 3 and 4 plot the correction to 
the energy density (p e g — p FRW ) / ' p FRW and the effective equation of state parameter w c s as a function of Ld at the 
current epoch. We have considered a set of Ld values equal to the simulation box size L divided by the powers of two, 
i.e. Ld = L, L/2, L/A, . . .. We can see from Fig. 3 that when Ld is comparable to the box size L, i.e. Ld = L, L/2, the 



t Note that is not the size of the whole domain, it is the smoothing length as described in [14-16]. 
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correction to the energy density is essentially negligible and we recover a standard FRW cosmology. This is because 
the local inhomogeneities in the individual zones are effectively smoothed out just as in the case when one simply 
uses the averaged energy density in a FRW cosmology. When Ld further decreases from L/4, one can clearly see that 
the correction quickly grows and then converges to a negative value at about the 10 -5 level. 

This can be explained by examining Eq. (33) for (p e g — Pfrw)/ Pfrw- It is obvious that the second term in this 
expression is always positive. Since our simulations are in the weak-held, slow-motion regime, the magnitude of the 
second term is much smaller than that of the first term because it is second order, i.e. <£> 2 . The sign of the first term 
is determined by the nature of the dominant volume weighted regions in the domain, i.e. collapsing overdense regions 
or expanding undcrdense regions. From Eq. (37), one can see that collapsing or expanding regions have negative or 
positive $ terms, respectively. This leads to a positive or negative first term. From Fig. 3, the simulation results 
clearly indicate that the expanding undcrdense regions are the dominant volume weighted regions. This also explains 
why the correction remains invariant as Ld further decreases. This is because once the resolution of the grid is fine 
enough to resolve the dominant regions, better resolution only slightly improves the accuracy of the correction. From 
Fig. 3 we can also find that the correction is nearly independent of the number of particles in the simulations. This 
is due to the way the matter density is distributed on the grid as described in Sec. III. 

The magnitude of the correction grows with the size of the simulation box. This is because simulations with larger 
box sizes have less restriction on the nonlinear growth of structures from a FRW cosmology boundary and thus include 
larger structures. It is clearly shown in Fig. 3 that the magnitude of the correction asymptotically converges to a 
value that is roughly represented by simulations with box sizes of 800 and 1600 Mpc/h. The most realistic asymptotic 
value is —5.6 x 10~ 5 at the best resolution of the grid in the simulations. 

For the equation of state parameter w e g, we only plot the results with Ld < L/4 in Fig. 4 because both P c g and 
PcS — Pfrw are close to zero for Ld = L,L/2, so that large numerical errors are introduced into the calculated w e g 
on these scales. It is shown in Fig. 4 that w e g is always negative and has a magnitude of a few tenths for the plotted 
Ld values. It decreases as Ld decreases. For different simulations, w e g increases with the box size and is almost 
independent of the number of particles in the simulations. Again, it asymptotically converges to a value given in the 
800 and 1600 Mpc/h simulations. The most realistic value is about -0.24. Note that w c g is negative because p e g is 
less than p FRW , while P c g is always positive. A cosmic acceleration requires a negative P c g. Therefore, even though 
W e fi < 0, no cosmic acceleration results. The sign of w e g and P e g can be explained by analyzing the expression for 
P ff , Eq. (40). Because our simulations are in the weak- field, slow- motion regime, all of the second order terms, e.g. 
$ 2 , <I> 2 and $<£", are much smaller than the first order terms. Hence, the first and third terms in the expression of P g 
are the dominant terms. Since the expanding undcrdense regions are the dominant regions in all of our simulations, 
the first term is positive. The third term is also positive. This is because the collapsing overdense regions' densities 
are always weighted more than the underdense ones due to the (1 — 2$) 3 factor no matter whether they are the 
dominant regions or not as long as 1 — 4$ 2 > 0. Therefore the effective pressure P e g is always positive. However, we 
will discuss one possible scenario in which one can have a large negative pressure term in Sec. V. 
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FIG. 3: Correction to the energy density, (p c g — Pfrw) / Pfrw as a function of the smoothing length Ld for various box sizes 
and numbers of particles. Note the convergence for the largest box sizes. 



As an illustration of the time evolution of the correction terms, Figs. 5 and 6 show (p c g — p FRW )/ Pfrw and 
w e g , respectively, as a function of redshift z with Ld = 3.125 Mpc/h. From Fig. 5, we can see that for all of the 
simulations, the correction to the energy density (p e g — p FRW ) / Pfrw grows from a negligible value at a redshift z ~ 40 
to a negative 10 -5 level at the current epoch. The raw data show that it grows by about five orders of magnitude 
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FIG. 4: Effective EOS parameter w e s as a function of the smoothing length Ld for various box sizes and numbers of particles. 
Note the convergence for the largest box sizes. 



during this process. Starting from a redshift z ~ 2, its magnitude quickly increases to the current level. The fact that 
it grows significantly during the structure formation process suggests it is the structure formation that leads to the 
correction to the energy density. In fact, this happens just as the dark energy becomes the dominant energy form in 
the Universe. This coincidence allows the possibility that there may yet be a possible connection between structure 
formation and the emergence of the dark energy. As before, the results represented by the simulations with box sizes 
of 800 and 1600 Mpc/h are the most realistic values. 

Fig. 6 shows that w c s is always negative and its magnitude grows gradually from a redshift z ~ 40 to the current 
epoch in all of the simulations. The apparent deviations, especially for the L = 1600 Mpc/h simulation, at high 
redshifts are due to the fact that both P c g and p c s — p FRW are very small at those redshifts and our simulations only 
have a limited number of particles and limited resolution, thus, some numerical errors are introduced into the values 
of w e g. This is not the case at low redshifts. For this reason and the trend we discussed before, we believe the result 
from the 800 Mpc/h with 512 3 particles simulation is closest to the true value of w e s- 
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FIG. 5: (p e s - Pfrw)/pf 
for the largest box sizes. 



as a function of redshift z for various box sizes and numbers of particles. Note the convergence 



From the results we have presented in this section, we conclude that: 1) Numerical simulations can be used to 
calculate the proposed effect despite previous claims to the contrary [31, 55-57]; 2) The correction to the energy 
density is negative and its magnitude is slightly larger than the previously claimed [37, 38] 10 -5 level at the current 
epoch. It grows by about five orders of magnitude along with the structure formation process; 3) The effective 
pressure is always positive and its magnitude is roughly of the same order as that of p e ff — Pfrw- Hence, the effective 
equation of state parameter is negative and its value is about -0.2 to -0.3; and 4) The proposed correction terms in the 
weak-field, slow-motion limit are not able to drive the current cosmic acceleration for the simple reasons that their 
magnitudes are too small and the effective pressure is always positive. 
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FIG. 6: w e fi as a function of redshift z for various box sizes and numbers of particles. Note the convergence for the largest box 
sizes. 



We have used the ADM formalism to develop a practical scheme to calculate the proposed domain averaging effect 
in an inhomogeneous cosmology within the context of numerical large-scale structure simulations. We find that in 
the weak-field, slow-motion limit, the proposed effect implies a small correction to the global expansion rate of the 
Universe. Under this limit, our simulations are always dominated by the expanding underdensc regions, hence the 
correction to the energy density is negative and the effective pressure is positive. However, whether this is still the case 
when strong-field gravity is included in a more general scenario needs to be further investigated. At least in the current 
investigation, the proposed effect cannot be the source of the current cosmic acceleration. We have done our analysis 
on each given redshift in standard FRW cosmology simulations, whether the cumulative effect can significantly change 
the expansion history of the Universe remains to be further studied. Nevertheless, the fact that this effect just begins 
to grow during the structure forming era allows the possibility that relativistic corrections from the development of 
the cosmic structure may have played a non-negligible role on the global dynamics of the Universe. 

We wish to point out that one possible scenario exists in which one could have a large negative correction to 
the energy density and a large negative pressure. This occurs when the underdense regions are the overwhelmingly 
dominant regions and they expand very rapidly during the cosmic evolution. In this scenario, the second order terms 
in the expression for the effective pressure become the dominant terms and they can induce a large negative pressure. 
We have verified this prediction in simple toy models. Whether this can happen in a realistic cosmological model 
needs to be investigated with strong-field gravity and full GR. If this is the case, the proposed GR correction in an 
inhomogeneous cosmology model may yet be found to serve as one possible source of the current cosmic acceleration. 
Also, for a simple collapsing overdense region in a FRW cosmology box, we find that both the correction to the energy 
density and the effective pressure are positive. In this scenario, the cosmic expansion is effectively slowed. From 
these two cases, it is suggested that the effect of the proposed GR correction on the local expansion rate behaves like 
a positive "feedback" on the structure formation. This aspect needs to be studied with a more realistic model and 
probably large-scale structure surveys. Unlike the conformal Newtonian gauge, the conformally flat model that we 
utilized in this paper can be applied to models beyond the weak-field, slow-motion limit. In future work, we will use 
this metric to investigate some of these aspects. 
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